perm filename SSM.F4[TMP,LCS] blob sn#340404 filedate 1975-04-18 generic text, type T, neo UTF8
	SUBROUTINE SSS(VV,N1,A1)
	DIMENSION V(50,4),A1(512),C(30,4),YP(30),J(30),NX(3),KA(14),K(9)
	DIMENSION VV(50,4)
	EQUIVALENCE(K1,K(1)),(K2,K(2)),(K3,K(3)),(K4,K(4)),(K5,K(5)),
     1	(K6,K(6)),(K7,K(7)),(K8,K(8)),(K9,K(9))
IDATA KA/1,2,2,1,1,2,1,1,0,2,1,-1,0,1/,DX/.00001/
	IF(VV(1,2).EQ.0) VV(1,2)=1
	DO 5 I=1,30
	DO 5 L=1,2
5	V(I,L)=VV(I,L)
	NX(1)=N1
698	NX(2)=NX(1)-1
	DO 10 I=1,NX(1)
10	V(I,2)=(V(I,2)-1)/99.
	DO 20 I=2,NX(2)
	JX=I+1
	JZ=I-1
	YP(I)=(V(JX,1)-V(JZ,1))/(V(JX,2)-V(JZ,2))
20	IF((V(JX,1)-V(I,1))*(V(I,1)-V(JZ,1)).LE.0) YP(I)=0
	DO 22 I=1,9
22	K(I)=KA(I)
	KOUNT=0
21	KOUNT=KOUNT+1
	V1=V(K2,1)-V(K1,1)
	V2=V(K2,2)-V(K1,2)
802	IF((YP(K2)-V1/V2)*(V(K3,1)-V(K4,1)).GT.0) GO TO 30
24	Z=V(K2,K5)+(V(K1,K6)-V(K2,K6))*YP(K2)**K7
	IF(YP(K2)**2.LT.DX.AND.V1**2.LT.DX) GO TO 36
	IF(YP(K2)**2.LT.DX) GO TO 38
	D1=V(K2,K5)-Z
806	D2=Z-V(K1,K5)
	ZZ=(V(K1,K6)*D2+V(K2,K6)*D1)/(D1+D2)
808	YP(K1)=(ZZ*K9+V(K2,1)*K8-V(K1,1))/
     1	(ZZ*K8+V(K2,2)*K9-V(K1,2))
	GO TO 40
30	DO 32 I=5,9
32	K(I)=KA(I+5)
	GO TO 24
36	YP(K1)=0
	GO TO 40
38	YP(K1)=-100
	IF(KOUNT.EQ.2) GO TO 39
	IF(V(K2,1).GT.V(K1,1)) YP(K1)=100
	GO TO 40
39	IF(V(K2,1).LT.V(K1,1)) YP(K1)=100
40	IF(KOUNT.EQ.2) GO TO 50
	DO 42 I=1,2
	K(I)=NX(I)
42	K(I+2)=K(I)     
	DO 44 I=5,9
44	K(I)=KA(I)
	GO TO 21
50	NX(3)=NX(2)-1
	N=1
52	N=N+1
	IF(N.GT.NX(3)) GO TO 92
	JX=N+1
	V1=V(JX,1)-V(N,1)
	V2=V(JX,2)-V(N,2)
	Y1=YP(N)-YP(JX)
	IF(Y1**2.LT.DX.AND.V1**2.GT.DX) GO TO 720
710	X=(V1-YP(JX)*V(JX,2)+YP(N)*V(N,2))/Y1                   
715	IF(X.GE.V(N,2).AND.X.LE.V(JX,2)) GO TO 52      
	IF(Y1**2.LT.DX.AND.V1**2.LT.DX) GO TO 52
720	DO 120 I=NX(1),JX,-1
	JZ=I+1
	V(JZ,2)=V(I,2)
	V(JZ,1)=V(I,1)
120	YP(JZ)=YP(I)
	YP(JX)=.5*V1/V2
	IF(V1*(YP(N)-V1/V2).LE.0) YP(N+1)=4*YP(JX)
	V(JX,2)=.5*(V(N+2,2)+V(N,2))
	V(JX,1)=.5*(V(N+2,1)+V(N,1))
	N=JX
	DO 88 L=1,3
88	NX(L)=NX(L)+1
  	GO TO 52
92	DO 140 I=1,NX(2)
	JX=I+1
	W0=YP(I)
	W1=YP(JX)
	W2=V(JX,2)-V(I,2)
	W3=V(JX,1)-V(I,1)
	C(I,1)=(W2*(W0+W1)-2*W3)/(W0-W1)
	C(I,2)=W2-C(I,1)
	C(I,4)=W0*C(I,2)
140	C(I,3)=-C(I,4)+W3
730	DO 150 I=1,NX(1)
150	J(I)=511*V(I,2)+1
740	DO 160 I=1,NX(2)
	L1=J(I)+1
	IF(I.EQ.1) L1=1
	ZZ=C(I,2)
	XX=C(I,1)
	L2=J(I+1)
750	DO 160 L=L1,L2
	X=(FLOAT(L)-1.)/511.
	IF(XX**2.LT.DX) GO TO 155
	ZX=.5*SQRT(ZZ**2-4*XX*(V(I,2)-X))/XX
	T1=-.5*ZZ/XX+ZX
	T2=T1-2*ZX
	IF(T2.GT.-DX.AND.T2.LT.(1+DX)) T1=T2
155	IF(XX**2.LT.DX) T1=-(V(I,2)-X)/ZZ
160   	A1(L)=C(I,3)*T1**2+C(I,4)*T1+V(I,1)
770	END